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ABSTRACT 

A' 1 ' 


The three dimensional viscous flow through a planar turbine cascade is numerically simulated by direct solution of the 
incompressible Navier-Stokes equations. Flow dependence in the spanwise direction is represented by direct expansion 
in Chebyshev polynomials, while the discretization on planes parallel to the endwalls is accomplished using the spectral 
element method. Elemental mapping from the physical to the computational space use an algebraic mapping technique. 
A fractional time stepping method that consists of an explicit nonlinear convective step, an implicit pressure correction 
step, and an implicit viscous step is used to advance the Navier-Stokes equations forward in time. Results computed 
at moderate Reynolds numbers show a three dimensional endwall flow separation, a midspan separation of the blade 
suction surface boundary layer, and other three-dimensional features such as the presence of a saddle point flow in the 
endwall region. In addition, the computed skin friction lines are shown to be orthogonal to the surface vorticity lines, 
demonstrating the accuracy achievable in the present method. 


INTRODUCTION 

Blade rows in modern axial flow turbines are often designed with rather high loadings and low aspect ratios (ref. [l]) 
(to increase power density and decrease part counts), this results in strong endwall secondary flows, often extending 
to midpan. These strong secondary flows have been known to be a source of (total pressure) loss through the blade 
row (ref. [2]). A large number of experimental investigations (ref. [3]-[4]) have been carried out in plane and annular 
cascades to obtain data which have been used to develop a physical understanding of the generation of secondary 
flows. Concurrent with these studies, theoretical efforts have been made (ref. [5]) to develop analytic models for 
predicting flowfields and the losses associated with the existence of secondary flows. These efforts have led to a global 
understanding of the secondary flow phenomena. However, the detailed mechanisms responsible for these secondary 
losses are not understood to the extent that a quantitative model can be formulated in terms of the secondary flow. 
Such an understanding is essential for developing a reliable model for the prediction of the losses associated with 
endwall secondary flows (ref. [6]). An approach to address this issue would be to make use of a numerical simulation 
scheme that can generate an accurate solution of the three-dimensional flow in a planar turbine cascade. The high 
order accuracy obtainable in a calculation based on the spectral element technique makes it ideally suited to such an 
investigation. 

This paper presents a computation of the three dimensional viscous flow through a planar turbine cascade by a multi- 
domain spectral element method. The goal of this investigation is to develop a reliable computational tool that can be 
used to gain an improved understanding of secondary flow and loss generation in turbine cascade endwall regions. The 
spectral element technique used in this investigation offers the advantages of high order accuracy, minimal dispersion 
and dissipation errors and geometric flexibility, all of which are essential to a quantitative study of the fundamental 
phenomena underlying the generation of secondary flow and the associated loss in turbine cascades. In the following, 
we first present the governing equations and outline an efficient technique for advancing the solution in time. This is 
followed by a brief concise description of the spatial discretization and formulation of the spectral element method. 
The last two sections present an application of the method to the calculation of the three dimensional viscous flow 
through a planar turbine cascade. 
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GOVERNING EQUATIONS AND TEMPORAL DISCRETIZATION 


The equations governing the flow are the incompressible Navier-Stokes equations written in rotational form, 


dV 

dt 


= V X w - + 



( 1 ) 


V-V=0 (2) 

Here, V is the velocity field normalized by upstream axial velocity at midspan, w = V x V is the vorticity field, P t is 
the total pressure normalized by twice the upstream axial dynamic head, and Re is the Reynolds number based on the 
upstream axial velocity at midspan and blade axial chord. 


The solution to Eqs.(l-2) is advanced forward in time using a fractional time stepping scheme (ref. [7]), consisting 
of a non-linear convective step, a pressure step imposing continuity, and a viscous correction step imposing the no-slip 
boundary condition. The non-linear convection step is implemented through an explicit third order Adams-Bashforth 
scheme that yields 

»n.+ l _ A + _ 

V -v n = — [23(V x w)" - 16(V x w)"" 1 + 5(V x w) n “ 2 ] (3) 


Once V is determined, we 
the pressure correction step is 


are left with an unsteady Stokes problem which can be split in time as follows, 
discretized in time by a Backward Euler method, yielding 


First, 


s n + l 
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a n+ 1 

- v 


At 


= -VP« 


( 4 ) 


and 


* n+l 

V-V = 0 


subject to the boundary condition 


v -4 = 0 


( 5 ) 

( 6 ) 


on the blade surface and endwalls. Computationally, the above step is reformulated as a solution for P t by taking the 
divergence of Eqn.(4) and applying Eqn.(5) to yield 


V 2 P t = V ■ ( 


- n+l 

V 

At ■ 


subject to the boundary condition 


i.n+1 

ap = v _4 

dn 


At 


( 7 ) 


( 8 ) 


^ n+l 


on the solid walls. The velocity field V that satisfies continuity identically is then computed from Eqn.(4). 


Following the solution of the pressure step is the viscous correction step imposing the non-slip boundary condition 
on the solid surfaces. The step is discretized using the implicit Crank- Nicholson scheme, giving 


(V 2 - ~)(V n+1 + V n ) 


2 Re 
~At 


2 n+l 

(v + v Tl ) 


( 9 ) 


subject to the appropriate non-slip boundary conditions on the solid surfaces. At the inflow, the velocity is assumed 
known, while at the outflow a homogeneous Neumann boundary condition is imposed. 
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SPATIAL DISCRETIZATION AND ELEMENTAL MAPPING 


Spatial discretization in Z 


Since the geometry is invariant in the spanwise dimension (planar cascade), one can choose a direct spectral 
expansion for the flow variation in the Z direction. The need to account for different boundary conditions on the 
endwalls (Dirichlet for the velocity and Neumann for the pressure) leads us to define an eigenfunction expansion 

F,{Z) =^f lm T m (Z) (10) 

m 


satisfying the following Sturm-Liouville problem, 


cPF,(Z) 

dZ 2 


= A fF,(Z), 


( 11 ) 


subject to homogeneous Neumann boundary conditions for the pressure and homogeneous Dirichlet boundary conditions 
for the velocity. These functions are constructed separately for the viscous velocity step and pressure step in a 
preprocessing procedure via the tau method (ref. [8]). The Chebyshev polynomials T m (Z) are given as 


T m (Z) = cos(m cos 1 Z) 


(12) 


with the collocation points chosen at 

Zk = cos——. (13) 

1 / 

It should be noted that the choice of collocation points can be arbitrary, yet the distribution given in Eqn.(13) can be 
shown to give an error that satisfies the minimax criterion (ref. [9]). In addition, such a choice also results in good 
resolution of the viscous boundary layers near the endwall. 


Spatial discretization in the X-Y plane 


The complexity of the geometry prohibits a simple global spectral discretization of the flow variables in the (X, Y) 
plane. The region is instead divided into a number of subdomains, or spectral elements, following the techniqe developed 
by Patera (ref. [10]). In each i th subdomain or element we can expand the flow variables as 



N x Ny Nz ( t7 n 


(rtlFiiZ) 


(14) 


where Fi(Z) are the interpolants from the direct expansion in the spanwise direction given in Eqn.( 10) above, and 
h m (S) are high order local Lagrangian interpolants in terns of Chebyshev polynomials. These can be writen as 


with 


2 M 1 

h-m(S) — ^ . _ T n (S m )T n (S) 

M n = 0 

(15) 

, j 1 for m yt 0 or M 

m (2 for m = 0 or M 

(16) 


where the S m are the collocation points in the computational space. 


Elemental Mapping 

The mapping from the phyical coordinate space (X, Y, Z) to the local natural coordinate space (£, r/, f) 1 is given by 
an isoparametric tensor-product mapping (ref. [ll]), 

J K 

( = ( 17 ) 

y=o k = 0 
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Figure 1: Element collocation grid in physical and computational space. 

where we have chosen ( = Z, with the collocation points in the elemental computational space defined as 

• / 7TJ ivk 7T IV 

(£>*?> *)yfci = I cos-jcos — cos— J 


(18) 


where j, k, l — 0 — > J, K, L. 

To complete the definition of the mapping from the physical coordinate space ( X , Y, Z) to the elemental computa- 
tional space (£, 77, f)\ we need to define the collocation points in the physical space (X, y)‘- fc . This can be accomplished 
by several different methods. The first is through the use of an analytic conformal functional mapping, which can be 
used whenever the elements are rectangular in some suitable regular curvilinear coordinate system. A second and less 
restrictive method uses a Laplace equation to generate a linear functional variation in two dimensions over the element 
in physical space. The functional values are then used to map the points in the computational space to the physical 
element. The third method is an algebraic method using elements with two linear and two generally curved sides. Due 
to its relative simplicity, ease of implementation and geometric generality, the algebraic method is used here. 

In the algebraic method, we first define a general parametric function X(S) = fx{S) and ^(5) = fy{S) on each 
side of each element. Collocation points are then distributed along each element side in arc length according to the 
formula 

| y } = fx,v(S m ) ( 19 ) 

' y m 

where the collocation arc lengths S m are defined as in Eqn.(l2). Next, making use of the linearity of two opposing 
sides, the interior points are defined along straight lines connecting the points on opposing curved sides, distributed 
according to Eqn.(19) as illustrated in Figure (l). 

The final discretized equations are obtained by substituting Eqn.(14) into the relevant temporal discretizations, 
Eqns.(3-9). For a complete derivation of the final discretized equations, and a detailed description of the computational 
cyle, see Tan (ref. [12]). 
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Figure 2: Elemental grid and cascade geometry 


APPLICATION 


Objective 

The computational method above has been used to calculate the three-dimensional viscous flow through a planar 
cascade at moderate Reynolds number. Although the flow through turbine cascades generally occurs at Reynolds 
numbers that would be considered “large”, a decision has been made to restrict the investigation to a laminar flow 
regime. This has several implications for the application of the results. The first is that the loss levels that are 
calculated will not be applicable to realistic machine geometries. The second is that any phenomena that has its origin 
in small scale turbulent motion will not be reproduced in the simulation. Lastly boundary layer parameters such as 
momentum and displacement thickness will not reflect those of the cascade operated at large Reynolds numbers. Yet 
it is felt that these restrictions do not invalidate the results of the simulation. Since the generation of secondary flow 
is a kinematic process with characteristic length scales of an order much larger than those associated with small scale 
turbulent motion, a laminar simulation will produce secondary flows with structure and form similar to the higher 
Reynold number case. In addition, the relationship between these secondary flow structures and the mechanisms by 
which they are responsible for an increase in loss should not depend on the small scale turbulent structure of the flow. 
Therefore, conclusions drawn about this relationship, based on the simulation, should be applicable to the general high 
Reynolds number turbulent flow situation. 

Cascade Geometry and Inflow Conditions 

The cascade blade cross section is that used by Langston (ref. [13]) in his benchmark experimental investigation, 
invariant in the spanwise direction, with a solidity of 1.0 and a blade aspect ratio of 2.0, both based on axial chord. The 
cascade extends over the region Z — (—1, 1) in the spanwise direction. Figure ( 2) shows the spatial discretization in 
the X , Y plane, where the X — Y plane is divided into 125 spectral elements, with each element representing a seven by 
seven by thirty three Chebyshev expansion, as indicated in Eqn.(14). Figure ( 2) also shows a graphic representation 
of one blade and one endwall, with each level on the blade surface corresponding to one collocation point of the direct 
expansion in the spanwise direction. 

The inflow condition to the cascade consist of a circumferentially uniform flow at a angle of 45.5 degrees from axial. 
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Figure 3: Midspan static pressure distribution (a), and wall shear stress distribution (b). 
The spanwise variation consists of an eighth order polynomial profile 
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representing the incoming endwall boundary layer. The boundary layer parameters for this layer are summarized in 
Table (1), normalized both by axial chord, and by leading edge radius. 


Table 1: Boundary layer parameters of inlet velocity profile 


Parameter 

(chord) 

(Rl.b.) 

S 

.5623 

10.46 

6 * 

.1110 

2.066 

e 

.0517 

.9612 

Re e 

51.666 


NUMERICAL RESULTS 

Presented in this section are numerical results for the flow through a planar turbine cascade, using the inflow 
conditions given above and a Reynolds number of 1000. The simulation was conducted with a time step size of 
At — 3.0 X 10~ 4 , which was limitted by stability restrictions on the explicit convection step. The results presented 
show the solution after six thousand time steps, or a nodimensional time of r = 1.8. 


Computed results of the flow at midspan are shown in Figures ( 3-5). Figure ( 3. a) show the midspan surface static 
pressure coefficient, defined as 
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Figure 4: Midspan total pressure loss coefficient contours. 


The distribution shown has a smooth chordwise variation, with none of the oscillations that are suggestive of insufficient 
resolution in a spectral calculation. This computed static pressure distribution agrees with that measured by Langston 
(ref. [13]), with the exception of a higher minimum suction side pressure. This suggests the presence of a laminar 
boundary layer separation on the suction surface. Indeed, the plot of shear stress distribution on the blade shown in 
Figure ( 3.b), indicates that this laminar separation occurs at a point on the suctions surface where the value of the 
wall shear stress vanishes, at X — .741. Positive shear stress on the suction surface behind the separation is indicative 
of reverse flow in the separated region. Contours of total pressure loss coefficient in Figure ( 4), defined as 


Cp T = 


Pr a> - Pt 
P T „ - Poo ' 


( 22 ) 


show a rapid thickening of the suction side boundary layer which then separates in the reqion of adverse pressure 
gradient. The comparitively thin boundary layer on the pressure side evolves with negligable total pressure loss up to 
mid-chord. This may suggest that only a small fraction of the blade profile loss is generated on the pressure surface. 
Examination of the results in Figure (4) indicate a region of low total pressure near the trailing edge. This region of 
low total pressure outside the boundary layer can be attributed to the unsteadiness associated with the presence of 
vortex shedding and the separated flow region. Contours of the static pressure distribution, shown in Figure (5), clearly 
indicate the occurence of vortex shedding. The closed circular countours immediately downstream of the trailing edge 
are characteristic of shed vorticies. 


Numerical results of the endwall flow region are presented in Figures (6-11). Figure (6) shows a comparison of 
the velocity and vorticity vectors at a height of A Z = .005 above the endwall. At this spanwise location very near 
the endwall the velocity and vortitity vectors are orthogonal and proportional, demonstrating the accuracy attainable 
using the described computational method. 


Figure (7. a) is a plot of the two dimensional projection of the velocity vectors near the endwall, at a height of 
A Z = .019 above the endwall. The figure clearly shows the separation line that is a result of the formation of the 
horseshoe vortex at the blade leading edge. The line begins just behind the saddle point, shown in Figure (7,b), and 
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Figure 5: Midspan static pressure coefficient contours. 




Figure 7: Velocity vectors above cascade endwall (a) and an enlargement of the region near the saddle point (b). 
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Figure 8: Total velocity contours at Z = —.98. 

extends across the blade passage toward the suction surface of the adjacent blade. Behind this line the velocity is 
signifigantly higher, indicating that the incoming endwall boundary layer separates and convects toward the blade 
suction surface, as can be seen in Figure (10), which shows contours of total pressure loss coefficient (Eqn.(22)) on 
constant axial planes. Figure (8) shows a contour plot of the total velocity magnitude at the same spanwise location 
as the vectors in Figure (7. a). The separation line can clearly be seen extending across the passage, with a signifigant 
increase in total velocity magnitude occuring across the line. This increase in velocity appears to be nearly uniform 
over the length to the separation line. Thus, the acceleration of the upstream boundary layer on one side of the line 
and the higher velocity fluid on the other side of the line are of equal magnitude. This means that the pressure gradient 
in the streamwise direction must be continuous across the line. This can be seen in Figure (9) which is a contour plot 
of the static pressure on the passage enwall. The plot shows that the strength of the streamwise pressure gradient is 
continuous accross the separation line, as the total velocity distribution suggested. The Figure also shows the “trough” 
of static pressure on the endwall that has been observed by Langston (ref. [13]). This trough is characteristic of the 
formation of a horseshoe vortex structure. 


Figure (11) shows contours of streamwise vorticity, defined as 



(23) 


plotted at a height A Z = .019 above the endwall. The figure shows the saddle point ahead of the horseshoe vortex and 
the formation of the vortex at the blade leading edge. The highest level of streamwise vorticity is concentrated directly 
behind the separation fine, underneath and in front of the horseshoe votex core. This structure can be identified as a 
counter vortex, cause by the high shear levels beneath the larger horseshoe vortex. 


CONCLUDING REMARKS 

A multi-domain spectral method is used to calculate the three-dimensional viscous flow through a planar turbine 
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Figure 9: Static pressure contours on cascade endwall. 

cascade at moderate Reynolds number. The solution uses a series expansion in Chebyshev polynomials to represent 
the flow dependence in the spanwise direction, while discretization on planes parallel to the endwalls is accomplished 
using the spectral element method. The elemental mapping from the physical to computational space is accomplished 
by using an algebraic mapping procedure developed for the calculation. The computational method is then applied to 
a planar cascade of turbine blades using Langston’s profile. 

The midspan static pressure loading is found to be qualitatively similar to that measured by Langston, with the 
exception of the presence of a boundary layer speration on the suction side. 

The chosen spatial discretization is found to give an adequate resolution of the flow features in the endwall region. 
These features include the formation of a horseshoe vortex about the leading edge, and the resulting separation of 
the upstream boundary layer. The static pressure distribution on the endwall shows the presence of a static pressure 
trough that is characteristic of the formation of a horseshoe vortex system. In addition, there is evidence of a small 
intense counter vortex underneath the horseshoe vortex. The separation of the incoming boundary layer results in the 
presence of high velocity fluid near the wall behind the separation line. This causes an increase in endwall shear stress 
in that region, which may result in increased total pressure loss in the cascade. 

The presented application of a spectral element technique to the calculation of turbine cascade flows shows that 
the spectral method can be used as a computational tool to gain an improved understanding of secondary flow in axial 
flow turbines. 
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Figure 11: Streamwise component of vorticity above endwall. 
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